Evaluating RNA localization model

Ideas discussed at last meeting:

1. Fit steady state model on different cells to evaluate whether \(\nu\), the transport bias parameter, is spatially varying across the tissue.
2. Aggregate into two boxes, 1 and 2.
3. Try forward simulating with different \(\nu\).

Does \(\nu\) vary spatially across the tissue?

First look at fitting models for different numbers of cells. With only oocyte and nurse cell 2 we get:

source('run_model_comparison_stst.R')
out = run_model_comparison_stst(identifier='Ststv203',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(1,2,rep(0,14)))

## Loading required package: ggplot2

## Loading required package: StanHeaders

## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## Loading required package: tidyr

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:rstan':
## 
##     extract

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

With only oocyte and neighbouring nurse cells (1,2,3,5,9):

out = run_model_comparison_stst(identifier='Ststv204',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(1,2,3,5,9,rep(0,11)))

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

With only oocyte and nurse cells within two teps from the oocyte (1,2,3,4,5,7,8,9,10,11,13):

out = run_model_comparison_stst(identifier='Ststv205',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(1,2,3,4,5,7,8,9,10,11,13,rep(0,5)))

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

And with all the nurse cells in the tissue:

out = run_model_comparison_stst(identifier='Ststv206',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = seq_len(16))

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

Just with the cells further away from the oocyte:

out = run_model_comparison_stst(identifier='Ststv207',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(rep(0,11),6,12,14,15,16))

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

Everything except the neighbouring nurse cells:

out = run_model_comparison_stst(identifier='Ststv208',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(rep(0,5),4,6,7,8,10,11,12,13,14,15,16))

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

Now see if we can summarise the distribution of values for \(\nu\) that we get for fitting to each different choice of cells

library(dplyr)
library(ggplot2)
stst_df = data.frame(nu=numeric(),xi=numeric(),phi=numeric(),which_cells=integer())
for (id in seq(from=3,to=8)){
  print(paste('fits/alt_save_MCStstv20',id,'.Rsave',sep=''))
  load(paste('fits/alt_save_MCStstv20',id,'.Rsave',sep='')) #load in the results of fitting at stst for each choice of cells
  e$which_cells = rep(id,length(e$nu))
  stst_df = full_join(stst_df,as.data.frame(e))  
}

## [1] "fits/alt_save_MCStstv203.Rsave"

## Joining, by = c("nu", "xi", "phi", "which_cells")

## Warning: Column `nu` has different attributes on LHS and RHS of join

## Warning: Column `xi` has different attributes on LHS and RHS of join

## Warning: Column `phi` has different attributes on LHS and RHS of join

## [1] "fits/alt_save_MCStstv204.Rsave"

## Joining, by = c("nu", "xi", "phi", "which_cells")

## Warning: Column `nu` has different attributes on LHS and RHS of join

## Warning: Column `xi` has different attributes on LHS and RHS of join

## Warning: Column `phi` has different attributes on LHS and RHS of join

## [1] "fits/alt_save_MCStstv205.Rsave"

## Joining, by = c("nu", "xi", "phi", "which_cells")

## Warning: Column `nu` has different attributes on LHS and RHS of join

## Warning: Column `xi` has different attributes on LHS and RHS of join

## Warning: Column `phi` has different attributes on LHS and RHS of join

## [1] "fits/alt_save_MCStstv206.Rsave"

## Joining, by = c("nu", "xi", "phi", "which_cells")

## Warning: Column `nu` has different attributes on LHS and RHS of join

## Warning: Column `xi` has different attributes on LHS and RHS of join

## Warning: Column `phi` has different attributes on LHS and RHS of join

## [1] "fits/alt_save_MCStstv207.Rsave"

## Joining, by = c("nu", "xi", "phi", "which_cells")

## Warning: Column `nu` has different attributes on LHS and RHS of join

## Warning: Column `xi` has different attributes on LHS and RHS of join

## Warning: Column `phi` has different attributes on LHS and RHS of join

## [1] "fits/alt_save_MCStstv208.Rsave"

## Joining, by = c("nu", "xi", "phi", "which_cells")

## Warning: Column `nu` has different attributes on LHS and RHS of join

## Warning: Column `xi` has different attributes on LHS and RHS of join

## Warning: Column `phi` has different attributes on LHS and RHS of join

stst_df = stst_df %>% mutate(which_cells=factor(which_cells))
h <- ggplot(data = stst_df,aes(x=nu,color=which_cells)) + geom_density() + theme_bw()
h

This could certainly be interpreted as evidence of non-uniform bias in transport through ring canals spatially across the tissue.

---

** Two super-compartments

To reduce stochastic variation due to small sample sizes, we further consider agreggating into two super-compartments. I assume well mixed concentrations of mRNA in each super-compartment:

\[\begin{array} \frac{\text{d}x_1} = 7a - b*(1-\nu)x_1 + b\nu x_2, \\ \frac{\text{d}x_1} = 8a - b\nu x_2 + b(1-\nu) x_1. \end{array}\]

From the long time limit, we get \[ x_2 = \frac{8a + b(1-\nu) x_1}{b\nu}. \]

library(rstan)
library(mvtnorm)
library(dplyr)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#read in data on WT
N=6 # sample size
data = matrix(as.numeric(read.csv('data/exp_data.csv',sep=',',header=FALSE,stringsAsFactors = FALSE)),ncol=16,byrow=TRUE)
data = data[seq_len(N),]
#aggregate into two compartments
x1_obs = data[,seq(from=1,to=16,by=2)] %>% apply(.,1,sum)
x2_obs = data[,seq(from=2,to=16,by=2)] %>% apply(.,1,sum)

  estimates <- stan(file = 'super_compartment_model.stan',
                    data = list (
                      N = N,
                      x1_obs = x1_obs,
                      x2_obs = x2_obs
                      ),
                    seed = 42,
                    chains = 4,
                    warmup = 1000,
                    iter = 2000,
                    control = list(adapt_delta = 0.99)
  )
print(estimates, pars = c('a','b','theta','sigma'))

## Inference for Stan model: super_compartment_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## a      34.62    0.74 26.65   1.50  13.67  28.63  49.36 101.01  1288    1
## b      79.57    1.32 61.30   3.79  31.63  64.86 114.72 225.94  2145    1
## theta   0.44    0.00  0.12   0.27   0.35   0.41   0.49   0.74  1241    1
## sigma 247.92    1.13 44.43 173.60 215.69 244.27 274.91 348.34  1547    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Aug 16 13:30:06 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

pairs(estimates, pars = c('a','b','theta','sigma'))

The results from the two super-compartment model contradict what we have found using the full data suggesting that transport through ring canals should be less biased than previously thought. However, these result rely upon an assumption of well mixed concentrations of mRNA within each super-compartment. This assumption does not hold. I would suggest that the results here may be misleading due to the inaccuracy of this assumption.

---

TODO: try forward simulating with different models, but particularly with different \(\nu\) for single connection to oocyte. Plus make predictions for WT and OE from full fits.

Another possibility: assume the oocyte is an absorbing state. Currently we assume bidirectional transport through each ring canal. We could hypothesize that the connections directly to the oocyte are `special’ in some way such that transport only occurs in one direction here. This should decrease the concentration of mRNA in the neighbouring nurse cells and increase in the oocyte in the model.

Might be plausible for the full dynamic model.

But at steady state, the absorbing oocyte case results in qualitatively different behaviour in the sense that for this case the eigenvector corresponding to the largest eigenvalue is [1,0 … 0] irrespective of the transport bias \(\nu\). mRNA accumulates in the oocyte but not in the nurse cells, whereas assuming all the ring canals are bidirectional results in a distribution of mRNA across cells dependent on \(\nu\) with nonzero values in the nurse cells.

library(rstan)
mc <- stan_model('model_comparison5.stan')
expose_stan_functions(mc) #get hold of functions defined in the stan code
nu = 0.72 #pick a vlaue for the bias from ring canals
B = construct_matrix(0.72)
set_oocyte_absorbing = function(B){
  B[,1] = 0
  return(B)
}
B = construct_matrix(nu)
B_absorb = set_oocyte_absorbing(B)

print(eigen(B))

## eigen() decomposition
## $values
##  [1] -1.917756e+00 -1.663930e+00 -1.493222e+00 -1.374166e+00 -1.202139e+00
##  [6] -1.163084e+00 -1.101435e+00 -1.000000e+00 -7.461741e-01 -7.102357e-01
## [11] -6.698615e-01 -6.258343e-01 -5.080965e-01 -4.627503e-01 -3.613152e-01
## [16]  1.976128e-16
## 
## $vectors
##               [,1]         [,2]        [,3]         [,4]        [,5]
##  [1,]  0.666340488 -0.663418548 -0.53797776 -0.664653592  0.46128567
##  [2,] -0.666340488 -0.107998368  0.53797776 -0.108199422 -0.46128567
##  [3,] -0.181435015  0.663418548 -0.32163923 -0.108199422 -0.12560151
##  [4,]  0.181435015  0.107998368  0.32163923 -0.017613859  0.12560151
##  [5,] -0.116039538  0.180639412  0.24850174  0.664653592  0.43121093
##  [6,]  0.116039538  0.029406416 -0.24850174  0.108199422 -0.43121093
##  [7,]  0.031595912 -0.180639412  0.14857102  0.108199422 -0.11741259
##  [8,] -0.031595912 -0.029406416 -0.14857102  0.017613859  0.11741259
##  [9,] -0.088193947  0.115530699  0.11894222  0.180975697 -0.18875355
## [10,]  0.088193947  0.018807323 -0.11894222  0.029461160  0.18875355
## [11,]  0.024013954 -0.115530699  0.07111164  0.029461160  0.05139490
## [12,] -0.024013954 -0.018807323 -0.07111164  0.004796003 -0.05139490
## [13,]  0.015358492 -0.031457362 -0.05494158 -0.180975697 -0.17644726
## [14,] -0.015358492 -0.005120966  0.05494158 -0.029461160  0.17644726
## [15,] -0.004181898  0.031457362 -0.03284776 -0.029461160  0.04804407
## [16,]  0.004181898  0.005120966  0.03284776 -0.004796003 -0.04804407
##              [,6]        [,7]         [,8]        [,9]       [,10]
##  [1,]  0.54381084  0.49517497  0.679901590 -0.42402411  0.49193414
##  [2,]  0.08852735 -0.49517497  0.110681654  0.42402411  0.08008230
##  [3,] -0.54381084  0.29604884  0.110681654  0.11545572 -0.49193414
##  [4,] -0.08852735 -0.29604884  0.018017944 -0.11545572 -0.08008230
##  [5,]  0.32512663  0.17451114  0.110681654  0.07384147 -0.13394665
##  [6,]  0.05292759 -0.17451114  0.018017944 -0.07384147 -0.02180527
##  [7,] -0.32512663  0.10433448  0.018017944 -0.02010598  0.13394665
##  [8,] -0.05292759 -0.10433448  0.002933154  0.02010598  0.02180527
##  [9,] -0.25119615 -0.28713514 -0.679901590 -0.52152792  0.45986119
## [10,] -0.04089240  0.28713514 -0.110681654  0.52152792  0.07486112
## [11,]  0.25119615 -0.17166866 -0.110681654  0.14200462 -0.45986119
## [12,]  0.04089240  0.17166866 -0.018017944 -0.14200462 -0.07486112
## [13,] -0.15018192 -0.10119308 -0.110681654  0.09082122 -0.12521365
## [14,] -0.02444822  0.10119308 -0.018017944 -0.09082122 -0.02038362
## [15,]  0.15018192 -0.06050002 -0.018017944 -0.02472932  0.12521365
## [16,]  0.02444822  0.06050002 -0.002933154  0.02472932  0.02038362
##             [,11]        [,12]       [,13]       [,14]       [,15]
##  [1,]  0.44367490  0.591241031 -0.46310664  0.56496562  0.55109740
##  [2,] -0.44367490  0.096248540  0.46310664  0.09197115 -0.55109740
##  [3,]  0.26525864  0.096248540  0.12609734 -0.56496562  0.32948302
##  [4,] -0.26525864  0.015668367 -0.12609734 -0.09197115 -0.32948302
##  [5,] -0.20494152 -0.591241031 -0.43291317  0.33777438  0.19421950
##  [6,]  0.20494152 -0.096248540  0.43291317  0.05498653 -0.19421950
##  [7,] -0.12252780 -0.096248540  0.11787609 -0.33777438  0.11611746
##  [8,]  0.12252780 -0.015668367 -0.11787609 -0.05498653 -0.11611746
##  [9,]  0.32668014  0.353483583 -0.18424064  0.19910698  0.15471424
## [10,] -0.32668014  0.057543839  0.18424064  0.03241276 -0.15471424
## [11,]  0.19531132  0.057543839  0.05016610 -0.19910698  0.09249856
## [12,] -0.19531132  0.009367602 -0.05016610 -0.03241276 -0.09249856
## [13,] -0.15089951 -0.353483583 -0.17222858  0.11903952  0.05452489
## [14,]  0.15089951 -0.057543839  0.17222858  0.01937853 -0.05452489
## [15,] -0.09021786 -0.057543839  0.04689539 -0.11903952  0.03259864
## [16,]  0.09021786 -0.009367602 -0.04689539 -0.01937853 -0.03259864
##               [,16]
##  [1,] -0.9490332012
##  [2,] -0.1544937769
##  [3,] -0.1544937769
##  [4,] -0.0251501497
##  [5,] -0.1544937769
##  [6,] -0.0251501497
##  [7,] -0.0251501497
##  [8,] -0.0040942104
##  [9,] -0.1544937769
## [10,] -0.0251501497
## [11,] -0.0251501497
## [12,] -0.0040942104
## [13,] -0.0251501497
## [14,] -0.0040942104
## [15,] -0.0040942104
## [16,] -0.0006664994

print(eigen(B_absorb))

## eigen() decomposition
## $values
##  [1] -1.8506469 -1.5891102 -1.4773990 -1.2839774 -1.1946481 -1.1423873
##  [7] -1.0921857 -0.8600000 -0.7384633 -0.6948672 -0.6649883 -0.5760226
## [13] -0.5002191 -0.4336353 -0.3414496  0.0000000
## 
## $vectors
##                [,1]          [,2]          [,3]          [,4]
##  [1,]  3.943199e-01 -4.505475e-01  4.150432e-01  5.366668e-01
##  [2,] -8.485429e-01 -4.571711e-15 -7.130051e-01  1.305241e-15
##  [3,]  3.227778e-16  8.325228e-01 -3.233952e-16 -1.398704e-16
##  [4,]  2.801942e-01 -1.001964e-15 -3.515094e-01 -3.824201e-16
##  [5,] -6.917153e-18  6.434661e-17  1.423595e-16 -8.012420e-01
##  [6,]  1.629329e-01  1.387343e-15  3.534887e-01  7.874473e-16
##  [7,] -5.719377e-17 -2.749042e-01 -4.177991e-17 -7.839515e-16
##  [8,] -5.380145e-02  1.134009e-15  1.742688e-01 -5.957913e-16
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  1.199176e-01  1.340689e-15  1.616794e-01 -4.123511e-16
## [11,] -4.485057e-17 -1.598567e-01  6.855239e-17  9.528992e-16
## [12,] -3.959755e-02  3.382269e-17  7.970747e-02  3.808785e-16
## [13,]  8.247258e-19 -1.243206e-17 -2.817746e-17  2.645751e-01
## [14,] -2.302596e-02  2.603537e-16 -8.015629e-02 -4.154894e-16
## [15,]  7.879598e-18  5.278570e-02  1.595270e-17  4.849430e-16
## [16,]  7.603318e-03 -9.846383e-16 -3.951681e-02  2.591564e-16
##                [,5]          [,6]          [,7]       [,8]          [,9]
##  [1,]  4.334712e-01  5.176061e-01 -4.982278e-01 -0.7071068 -5.799736e-01
##  [2,] -6.021460e-01 -1.377579e-14  6.327410e-01  0.0000000  4.980107e-01
##  [3,]  2.960621e-15 -6.875658e-01 -3.279721e-16  0.0000000  2.987860e-16
##  [4,]  1.988324e-01 -7.902666e-16  3.119394e-01  0.0000000 -1.644462e-01
##  [5,]  1.782552e-16  1.590424e-16  3.467737e-18  0.0000000 -9.722521e-18
##  [6,] -5.105009e-01 -7.471665e-15  2.077652e-01  0.0000000 -9.562545e-02
##  [7,]  1.187878e-15 -3.389679e-01 -2.080447e-16  0.0000000 -3.342942e-17
##  [8,]  1.685706e-01  1.748434e-15  1.024276e-01  0.0000000  3.157613e-02
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.7071068  0.000000e+00
## [10,]  2.519077e-01  5.936144e-15 -3.815211e-01  0.0000000  5.736660e-01
## [11,] -1.240031e-15  3.408765e-01  2.133050e-16  0.0000000  2.289866e-16
## [12,] -8.318151e-02  4.790898e-16 -1.880887e-01  0.0000000 -1.894281e-01
## [13,] -9.758457e-17 -9.307494e-18 -4.009275e-18  0.0000000  1.332493e-17
## [14,]  2.135680e-01  3.309307e-15 -1.252753e-01  0.0000000 -1.101524e-01
## [15,] -5.446099e-16  1.680511e-01  1.113667e-16  0.0000000  8.644795e-17
## [16,] -7.052150e-02 -1.089499e-16 -6.176032e-02  0.0000000  3.637301e-02
##               [,10]         [,11]         [,12]         [,13]
##  [1,]  6.674936e-01  6.450812e-01 -8.012420e-01  7.575452e-01
##  [2,] -8.803343e-15 -4.988041e-01  1.398727e-16 -4.406262e-01
##  [3,] -5.393249e-01 -3.306966e-17  6.086083e-16  1.434219e-15
##  [4,]  1.320414e-15 -2.459089e-01  6.343025e-16  1.454976e-01
##  [5,] -6.896561e-17 -1.493518e-16  5.366668e-01 -3.212961e-16
##  [6,] -7.532866e-16  2.472936e-01 -5.930919e-16 -3.735641e-01
##  [7,]  1.780885e-01 -2.218286e-16 -5.989054e-16  4.644482e-16
##  [8,]  2.087871e-15  1.219150e-01 -8.141743e-17  1.233532e-01
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,] -6.653132e-15 -3.580943e-01  7.763897e-17 -1.714590e-01
## [11,] -4.572410e-01  7.301365e-17  9.532322e-16  5.286705e-16
## [12,]  5.623556e-17 -1.765394e-01  4.809293e-16  5.661684e-02
## [13,]  4.604350e-17 -1.450831e-16  2.645751e-01 -1.067148e-16
## [14,]  1.040675e-15  1.775335e-01 -3.130082e-16 -1.453634e-01
## [15,]  1.509839e-01 -2.285894e-16 -3.722299e-16  2.178551e-16
## [16,]  7.177216e-16  8.752347e-02 -1.399036e-16  4.799990e-02
##               [,14]         [,15] [,16]
##  [1,]  8.606270e-01 -9.006089e-01     1
##  [2,]  6.511794e-15  3.575727e-01     0
##  [3,] -4.339514e-01 -3.215760e-16     0
##  [4,] -3.316949e-16  1.762823e-01     0
##  [5,] -2.695983e-16 -1.068387e-17     0
##  [6,]  4.292391e-15  1.174116e-01     0
##  [7,] -2.139367e-01 -1.348413e-16     0
##  [8,] -9.261474e-16  5.788358e-02     0
##  [9,]  0.000000e+00  0.000000e+00     0
## [10,]  2.344028e-15  9.653869e-02     0
## [11,] -1.424911e-01 -8.822899e-17     0
## [12,] -4.093421e-16  4.759328e-02     0
## [13,] -6.069709e-17 -1.802243e-18     0
## [14,]  1.762836e-15  3.169919e-02     0
## [15,] -7.024770e-02 -3.379160e-17     0
## [16,] -5.409596e-16  1.562760e-02     0

Now try to visualise the effect of changing \(\nu\) by forward simulating from the model. (For this we will assume bidirectional transport through all ring canals again).

library(tidyr)
#keep all other parameters constant during this experiment
m0 = c(0, rep(1,15)) #initial condition
th = c(4,300) #c(96.5,249)
sig = 50
phi = 0.289
nSamples = 15
nTest = 5
nTotal = nSamples + nTest
source('extract_times_and_scaling.R')
times = extract_times_and_scaling(nSamples,nTest,test_on_mutant_data=FALSE)
data = matrix(as.numeric(read.csv('data/exp_data.csv',sep=',',header=FALSE,stringsAsFactors = FALSE)),ncol=16,byrow=TRUE)
raw_data = data[times$sort_indices2,]
#take a vector of transport bias values for nu to loop through
nu_vec = seq(from=0.7, to=1, by=0.05)
all_extracted_samples = data.frame(rna=numeric(),time=numeric())
for (j in seq_along(nu_vec)){
B = construct_matrix(nu_vec[j])
samples <- stan(file = 'mrna_transport6.stan',
                  data = list (
                    T  = nTotal,
                    y0 = m0,
                    t0 = times$t0,
                    ts = times$ts2,
                    theta = array(th, dim = 2),
                    sigma = sig,
                    phi = phi,
                    B = B %>% as.vector()
                  ),
                  algorithm="Fixed_param",
                  seed = 42,
                  chains = 1,
                  iter =100, 
                  refresh = -1
  )
estimates = rstan::extract(samples, pars='y_hat')
#make the samples go in a 'nice' format
pred <- as.data.frame(estimates, pars = 'y_hat') %>%
  gather(factor_key = TRUE) %>%
  group_by(key) %>%
  summarize(lb = quantile(value, probs = 0.05),
            median = quantile(value, probs = 0.5),
            ub = quantile(value, probs = 0.95))
xdata <- data.frame(rna = as.vector(raw_data),cellID = as.vector(matrix(rep(1:16,nTotal),nrow=nTotal,byrow=TRUE)),time = rep(times$ts2,16))
pred <- pred %>% bind_cols(xdata) %>%
        mutate(split = if_else(time %in% times$ts3,'train','test'),
               nu = nu_vec[j])
all_extracted_samples = full_join(all_extracted_samples, pred)
}

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.04661 seconds (Sampling)
##                0.04661 seconds (Total)

## Joining, by = c("rna", "time")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.04929 seconds (Sampling)
##                0.04929 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.052461 seconds (Sampling)
##                0.052461 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.056709 seconds (Sampling)
##                0.056709 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.054649 seconds (Sampling)
##                0.054649 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.059533 seconds (Sampling)
##                0.059533 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.055752 seconds (Sampling)
##                0.055752 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

p1 <- ggplot(all_extracted_samples, aes(x = time, y = median, group = factor(nu), color=factor(nu)))
p1 <- p1 + geom_line() +
  facet_wrap(~cellID,scales='free') +   #needed to remove factor(cellID) 
  labs(x = "time", y = "rna") +
  theme(text = element_text(size = 12), axis.text = element_text(size = 12),
        strip.text = element_text(size = 8)) + 
  scale_color_discrete(name = "nu") + 
  scale_colour_brewer(palette=7) + 
  geom_point(aes(x = time, y = rna))

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

p1

## Warning: Removed 1078 rows containing missing values (geom_point).

And finally compare this to what happens when we change one the transport bias for ring canals to the oocyte, with bias \(\nu\) for the other ring canals fixed at \(0.9\).

B1 = construct_matrix(0.72) #we fix ring canals for all other nurse cells
#take a vector of transport bias values for nu to loop through
nu_vec = seq(from=0.7, to=1, by=0.05)
all_extracted_samples = data.frame(rna=numeric(),time=numeric())
for (j in seq_along(nu_vec)){
#mess around with ring canals into oocyte 
  B2 = construct_matrix(nu_vec[j])
  B = B1
  B[,1] = B2[,1]
  B[1,] = B2[1,]
  diag(B) = diag(B) - apply(B,2,sum)
samples <- stan(file = 'mrna_transport6.stan',
                  data = list (
                    T  = nTotal,
                    y0 = m0,
                    t0 = times$t0,
                    ts = times$ts2,
                    theta = array(th, dim = 2),
                    sigma = sig,
                    phi = phi,
                    B = B %>% as.vector()
                  ),
                  algorithm="Fixed_param",
                  seed = 42,
                  chains = 1,
                  iter =100, 
                  refresh = -1
  )
estimates = rstan::extract(samples, pars='y_hat')
#make the samples go in a 'nice' format
pred <- as.data.frame(estimates, pars = 'y_hat') %>%
  gather(factor_key = TRUE) %>%
  group_by(key) %>%
  summarize(lb = quantile(value, probs = 0.05),
            median = quantile(value, probs = 0.5),
            ub = quantile(value, probs = 0.95))
xdata <- data.frame(rna = as.vector(raw_data),cellID = as.vector(matrix(rep(1:16,nTotal),nrow=nTotal,byrow=TRUE)),time = rep(times$ts2,16))
pred <- pred %>% bind_cols(xdata) %>%
        mutate(split = if_else(time %in% times$ts3,'train','test'),
               nu = nu_vec[j])
all_extracted_samples = full_join(all_extracted_samples, pred)
}

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.057958 seconds (Sampling)
##                0.057958 seconds (Total)

## Joining, by = c("rna", "time")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.052128 seconds (Sampling)
##                0.052128 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.0512 seconds (Sampling)
##                0.0512 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.051745 seconds (Sampling)
##                0.051745 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.052531 seconds (Sampling)
##                0.052531 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.049678 seconds (Sampling)
##                0.049678 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.049038 seconds (Sampling)
##                0.049038 seconds (Total)

## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")

p2 <- ggplot(all_extracted_samples, aes(x = time, y = median, group = factor(nu), color=factor(nu)))
p2 <- p2 + geom_line() +
  facet_wrap(~cellID,scales='free') +   #needed to remove factor(cellID) 
  labs(x = "time", y = "rna") +
  theme(text = element_text(size = 12), axis.text = element_text(size = 12),
        strip.text = element_text(size = 8)) + 
  scale_color_discrete(name = "nu") + 
  scale_colour_brewer(palette=7) + 
  geom_point(aes(x = time, y = rna))

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

p2

## Warning: Removed 1078 rows containing missing values (geom_point).

Changing the transport bias, \(\nu\), for selected ring canals neighbouring the oocyte has limited effect relative to predictions from the simpler model where \(\nu\) is assumed to be spatially uniform.

Seems possible though to look at fitting this model with spatially varying \(\nu\) values based on whether cells are neighbouring the oocyte. Can then compare between the models to assess whether this is plausible.

source('mrna_transport_full.R')
mrna_transport_inference_full('full_spatial_nu_v201', use_real_data = TRUE, run_mcmc = FALSE, nSamples=6, nTest=14, show_diagnostic_plots = TRUE, test_on_mutant_data = FALSE)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]  566   13   18   NA   16   NA   NA   NA   10    NA    NA    NA    NA
##  [2,]  708   47  117   50   91   34   86   49   76    94   100    63   110
##  [3,]  800   44  119   NA   29   NA   NA   NA    5    NA    NA    NA    NA
##  [4,] 1106  162   36  118   92  131   61   76   82   129    39   135   122
##  [5,]  446   83   22   NA    2   NA   NA   NA   44    NA    NA    NA    NA
##  [6,]  982   77   69   NA   15   NA   NA   NA    0    NA    NA    NA    NA
##  [7,] 1185   75  122   NA   51   NA   NA   NA   23    NA    NA    NA    NA
##  [8,] 1076   22   19   NA   67   NA   NA   NA   70    NA    NA    NA    NA
##  [9,] 1007   51    5   NA   32   NA   NA   NA   11    NA    NA    NA    NA
## [10,] 1032   31  144   NA   43   NA   NA   NA    6    NA    NA    NA    NA
## [11,] 1078   47  159   NA   18   NA   NA   NA   61    NA    NA    NA    NA
## [12,] 1454  261  144  180  108  142   10   96   18   127    20    53    19
## [13,] 1012   23  141   NA   28   NA   NA   NA   88    NA    NA    NA    NA
## [14,] 1862   48   32   NA   73   NA   NA   NA   21    NA    NA    NA    NA
## [15,] 1045   78  130   NA  101   NA   NA   NA   11    NA    NA    NA    NA
## [16,] 1840  283  194  168  190  179   67   22   85    70    56    47   174
## [17,] 1365  102  318   NA  155   NA   NA   NA   31    NA    NA    NA    NA
## [18,] 1615   98   80   NA  284   NA   NA   NA  164    NA    NA    NA    NA
## [19,] 1485  407   53  384  399  174   47  152   21   315    37   176   266
## [20,] 1652   68  376  249   43  166   28  330  263   293   208    73    29
##       [,14] [,15] [,16]
##  [1,]    NA    NA    NA
##  [2,]    73    81    71
##  [3,]    NA    NA    NA
##  [4,]    84    31    35
##  [5,]    NA    NA    NA
##  [6,]    NA    NA    NA
##  [7,]    NA    NA    NA
##  [8,]    NA    NA    NA
##  [9,]    NA    NA    NA
## [10,]    NA    NA    NA
## [11,]    NA    NA    NA
## [12,]   109    20    11
## [13,]    NA    NA    NA
## [14,]    NA    NA    NA
## [15,]    NA    NA    NA
## [16,]   136    20    17
## [17,]    NA    NA    NA
## [18,]    NA    NA    NA
## [19,]    40    99    73
## [20,]   272    12    98

## Warning: Removed 154 rows containing missing values (geom_point).

## Saving 7 x 5 in image

## Warning: Removed 154 rows containing missing values (geom_point).

## Joining, by = "parameter"

## Saving 7 x 5 in image

## This is bayesplot version 1.6.0

## - Online documentation and vignettes at mc-stan.org/bayesplot

## - bayesplot theme set to bayesplot::theme_default()

##    * Does _not_ affect other ggplot2 plots

##    * See ?bayesplot_theme_set for details on theme setting

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## the following parameters were dropped because they are duplicative
## a b

## Warning: Removed 154 rows containing missing values (geom_point).

Finally, make predictions for full model for WT and OE, with spatially uniform \(\nu\) and spatially varying \(\nu\).